/*
* Copyright 2003 Program of Computer Graphics, Cornell University
*     580 Rhodes Hall
*     Cornell University
*     Ithaca NY 14853
* Web: http://www.graphics.cornell.edu/
* 
* Not for commercial use. Do not redistribute without permission.
*/

// This file contains a simple testing harness for the Brite ray engine
// It can load simple geometry files (.rrf format) and ray files (.rrf or .rrfb)
// and intersect the loaded rays against the loaded geometry
//

//use slower but higher precision triangle intersection tests
//#define DOUBLE_PRECISION_TRIANGLES
#include "stdafx.h"

// Includes for the timer
#ifdef _WIN32
#define WIN32_LEAN_AND_MEAN		// Exclude rarely-used stuff from Windows headers
#include <windows.h>
#else
#include <sys/time.h>
#include <unistd.h>
#endif
#include <cstdlib>

#include "Ray.h"
#ifdef DOUBLE_PRECISION_TRIANGLES
#include "TriangleDP.h"
#else
#include "Triangle.h"
#endif
#include "Parallelogram.h"
#include "Scene.h"

class RayWithResult {
public:
	float originX,originY,originZ;
	float dirX,dirY,dirZ;
	float tStart,tEnd,tResult;
	void constructRay(Ray &ray) const {
		ray.setOrigin(originX,originY,originZ);
		ray.setDirection(dirX,dirY,dirZ);
		ray.setTParameters(tStart,tEnd);
		assert(originX==ray.originX());
	}
};

// Timer with sub-millisecond accuracy.
// Taken from the Portable Agile C++ Classes (PACC)
// http://manitou.gel.ulaval.ca/~parizeau/PACC
class Timer {

protected:
	bool mHardware; 
	unsigned long long mCount; 
	unsigned long long elapsed;
	static double mPeriod; 

public:
	Timer(bool inHardware=true) : mHardware(inHardware) {
		if(mPeriod == 0) calibrateCountPeriod(); 
		reset();
	}

	void calibrateCountPeriod(unsigned int inDelay=10000, unsigned int inTimes=10)
	{
#ifdef _WIN32
		// use the windows counter
		LARGE_INTEGER lFrequency;
		QueryPerformanceFrequency(&lFrequency);
		mPeriod = 1. / lFrequency.QuadPart;
#else
		if(mHardware) {
#if defined (__GNUG__) && (defined (__i386__) || defined (__ppc__))
			double lPeriod = 0;
			// calibrate by matching the time-stamps with the micro-seconds of gettimeofday
			for(unsigned int i = 0; i < inTimes; ++ i) {
				timeval lStartTime, lTime;
				::gettimeofday(&lStartTime, 0);
				unsigned long long lStartCount = getCount();
				::usleep(inDelay);
				::gettimeofday(&lTime, 0);
				unsigned long long lCount = getCount() - lStartCount;
				lTime.tv_sec -= lStartTime.tv_sec;
				lTime.tv_usec -= lStartTime.tv_usec;
				// dismiss the first run of the loop
				if(i != 0) lPeriod += (lTime.tv_sec + lTime.tv_usec*0.000001)/lCount;
			}
			mPeriod = lPeriod/(inTimes-1);
#else
			// use the microseconds of gettimeofday
			mPeriod = 0.000001;
#endif
		} else {
			// use the microseconds of gettimeofday
			mPeriod = 0.000001;
		}
#endif

	}

	unsigned long long getCount(void) const
	{
		unsigned long long lCount = 0;
#ifdef _WIN32
		LARGE_INTEGER lCurrent;
		QueryPerformanceCounter(&lCurrent);
		lCount = lCurrent.QuadPart;
#else
		if(mHardware) {
#if defined (__GNUG__) && defined (__i386__)
			__asm__ volatile("rdtsc" : "=A" (lCount));
#else
#if defined (__GNUG__) && defined (__ppc__)
			register unsigned int lLow;
			register unsigned int lHigh1;
			register unsigned int lHigh2;
			do {
				// make sure that high bits have not changed
				__asm__ volatile ( "mftbu %0" : "=r" (lHigh1) );
				__asm__ volatile ( "mftb  %0" : "=r" (lLow) );
				__asm__ volatile ( "mftbu %0" : "=r" (lHigh2) );
			} while(lHigh1 != lHigh2);
			// transfer to lCount
			unsigned int *lPtr = (unsigned int*) &lCount;
			*lPtr++ = lHigh1; *lPtr = lLow;
#else
			timeval lCurrent;
			::gettimeofday(&lCurrent, 0);
			lCount = (unsigned long long)lCurrent.tv_sec*1000000 + lCurrent.tv_usec;
#endif
#endif
		} else {
			timeval lCurrent;
			::gettimeofday(&lCurrent, 0);
			lCount = (unsigned long long)lCurrent.tv_sec*1000000 + lCurrent.tv_usec;
		}
#endif
		return lCount;

	}

	double getCountPeriod(void) const {return mPeriod;}

	void reset(void) {mCount = getCount(); elapsed=0;}
	void start(void) {mCount = getCount();}
	void stop(void)  {elapsed += getCount()-mCount;}
	double elapsedSeconds() {return elapsed*mPeriod;} 
};
double Timer::mPeriod = 0.0;


// fopen_s and sscanf_s seem to be available only in MSVC 2005,
// so we provide fixes for other compilers
#if !defined(_MSC_VER) || (_MSC_VER+0 < 1400)
inline int fopen_s(FILE **f, const char *filename, const char *mode) {
	*f = fopen(filename, mode);
	return f != NULL ? 0 : 1;
}
#define sscanf_s sscanf
#endif


bool loadGeometryFromFile(const char *filename, Scene *scene)
{
	int numTris = 0;
	int numOther = 0;
	int line = 0;
	int id = 0;
	int err;
	FILE *f;
	err = fopen_s(&f, filename, "r");
	if (err != 0) {
		printf("Unable to open file %s\n",filename);
		return false;
	}
	PoolAllocator *geomPool = scene->getGeometryObjectPool();
	char buf[4096];
	while(fgets(buf,sizeof(buf),f)) {
		line++;
		if (buf[0]=='t') {
			float v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z;
			if (sscanf_s(buf,"tri %g %g %g %g %g %g %g %g %g",
				&v0x,&v0y,&v0z,&v1x,&v1y,&v1z,&v2x,&v2y,&v2z) != 9) {
					printf("bad triangle data on line %d: %s\n",line,buf);
					return false;
			}
#ifdef DOUBLE_PRECISION_TRIANGLES
			TriangleDP *tri = TriangleDP::poolAlloc(geomPool);
			tri.set(v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,id++);
#else
			Triangle *tri = Triangle::poolAlloc(geomPool);
			tri->set(v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,id++);
#endif
			scene->addGeometryObject(tri);
			numTris++;
		} else if (buf[0]=='%') {
			//its a comment so skip it
		} else if (buf[0]=='r') {
			float v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z;
			if (sscanf_s(buf,"rect %g %g %g %g %g %g %g %g %g %g %g %g",
				&v0x,&v0y,&v0z,&v1x,&v1y,&v1z,&v2x,&v2y,&v2z,&v3x,&v3y,&v3z) != 12) {
					printf("bad rectangle data on line %d: %s\n",line,buf);
					return false;
			}
			Parallelogram *para = Parallelogram::poolAlloc(geomPool);
			para->set(v0x,v0y,v0z,v1x,v1y,v1z,v2x,v2y,v2z,v3x,v3y,v3z,id++);
			scene->addGeometryObject(para);
			numOther++;
		} else {
			printf("unrecognized data on line %d: %s\n",line,buf);
			return false;
		}
	}
	printf("%d triangles and %d others loaded from %s\n",numTris,numOther,filename);
	return true;
}

static inline void endianByteSwapFloats(int count, float *data) {
#if !(defined (__GNUG__) && defined (__ppc__))
	char *p = (char *)data;
	for(int i=0;i<count;i++) {
		char a = p[0], b = p[1], c = p[2], d = p[3];
		p[0] = d;
		p[1] = c;
		p[2] = b;
		p[3] = a;
		p += 4;
	}
#endif
}

bool loadRaysFromFile(const char *filename, ExpandableArray<RayWithResult> &rays)
{
	const float extraEpsilon = 2.0e-3f;
	//  const float extraEpsilon = 5.0e-3;
	int numRays = 0;
	int line = 0;
	//int id = 0;
	int err;
	FILE *f;
	err = fopen_s(&f, filename, "rb");
	if (err != 0) {
		printf("Unable to open file %s\n",filename);
		return false;
	}
	//before we start parsing the file, compute its size
	fseek(f, 0, SEEK_END);
	long fileSize=ftell(f);
	rewind(f);
	rays.removeAll();
	//now look for binary magic string at the start of the file
	char magic[16];
	memset(magic,0,16);
	if (fread(magic,16,1,f) != 1) {
		printf("Unable to read magic number for file %s\n",filename);
		return false;
	}
	if (strncmp(magic,"RayFileBinaryPCG",16)!=0) {
		rewind(f); //its not binary, so read as ascii
		char buf[4096];
		while(fgets(buf,sizeof(buf),f)) {
			line++;
			if (buf[0]=='r') {
				RayWithResult r;
				if (sscanf_s(buf,"ray %g %g %g %g %g %g %g %g %g",&r.originX,&r.originY,&r.originZ,
					&r.dirX,&r.dirY,&r.dirZ,&r.tStart,&r.tEnd,&r.tResult)!=9) {
						printf("bad ray data on line %d: %s\n",line,buf);
						return false;
				}
				if (r.tEnd >= 1e37f) r.tEnd = FLT_MAX;
				if (r.tResult >= 1e37f) r.tResult = FLT_MAX;
				r.tStart += extraEpsilon;
				r.tEnd -= extraEpsilon;
				if (r.tEnd < r.tStart) r.tEnd = r.tStart;
				rays.append(r);
				numRays++;
			} else if (buf[0]=='%') {
				//its a comment so skip it
			} else {
				printf("unrecognized data on line %d: %s\n",line,buf);
				return false;
			}
		}
	} else {
		rays.ensureCapacity(int(fileSize)/sizeof(RayWithResult));
		for(;;) {
			RayWithResult r;
			if (fread(&r,sizeof(RayWithResult),1,f)!=1) {
				if (ferror(f)) {
					printf("bad ray data\n");
					return false;
				} else {
					break;
				}
			}
			endianByteSwapFloats(9,(float*)&r);
			if (r.tEnd >= 1e37) r.tEnd = FLT_MAX;
			if (r.tResult >= 1e37) r.tResult = FLT_MAX;
			r.tStart += extraEpsilon;
			r.tEnd -= extraEpsilon;
			if (r.tEnd < r.tStart) r.tEnd = r.tStart;
			rays.append(r);
			numRays++;
		}
	}
	printf("%d rays loaded from %s\n",numRays,filename);
	return true;
}

void shootRays(Scene *scene, const ExpandableArray<RayWithResult> &rays, float *results)
{
	Timer timer;
	Ray ray;
	IntersectionRecord hitRec;
	int count = rays.size();
	timer.start();
	for(int i=0;i<count;i++) {
		hitRec.primitiveId = -1;
		rays[i].constructRay(ray);
		//    scene->intersect(ray,&hitRec);
		scene->intersect(1,&ray,&hitRec);
		results[i] = (hitRec.primitiveId>=0)?hitRec.t:FLT_MAX;
	}
	timer.stop();
	double secs = timer.elapsedSeconds();
	printf("%g seconds for %d rays (%g us per ray)\n",secs,count,1e6*secs/count);
}

float analyzeResults(Scene *scene, const ExpandableArray<RayWithResult> &rays, float *results, bool verbose)
{
	int raysGood = 0;
	int count = rays.size();
	for(int i=0;i<count;i++) {
		float error = absv(rays[i].tResult-results[i]) / (absv(rays[i].tResult)+absv(results[i]));
		if (error > 2e-4) {
			if (verbose)
				printf("rel error on ray %d: %g (%g vs %g)\n",i,error,rays[i].tResult,results[i]);
		} else {
			raysGood++;
		}
	}
	float successRatio = 100.0f*raysGood/count;
	printf("%d of %d rays good (%g %%)\n",raysGood,count,successRatio);
	return successRatio;
}

void buildAccelerationStructure(Scene *scene) 
{
	Timer timer;
	timer.start();
	//    scene->buildFlatList();
	//    scene->buildNaiveMiddleSplit(15);
	//    scene->buildNaiveMedianSplit(8);
	scene->buildGreedySplit();
	timer.stop();
	double secs = timer.elapsedSeconds();
	printf("%g seconds to build acceleration structure\n",secs);
	scene->printTreeStats();
}

int main(int argc, char* argv[])
{
	// Get the scene and rays files
	if (argc != 3) {
		fprintf(stderr, 
			"Incorrect number of arguments. Usage:\n  %s <scene_file> <rays_file>\n", 
			argv[0]);
		return EXIT_FAILURE;
	}
	const char *sceneFile = argv[1];
	const char *raysFile  = argv[2];

	puts("BRITE ray tester...");
	Scene *scene = new Scene();
	ExpandableArray<RayWithResult> rays;
	//    if (!loadGeometryFromFile("C:/raydata/box_lmbrtn.rrf",scene)) return -1;
	//    if (!loadRaysFromFile("C:/raydata/boxrays128_1.rrfb",rays)) return -1;
	//    if (!loadGeometryFromFile("C:/raydata/taliesin.rrf",scene)) return -1;
	//    if (!loadRaysFromFile("C:/raydata/taliesinrays128_1.rrfb",rays)) return -1;
	//    if (!loadGeometryFromFile("C:/raydata/bar.rrf",scene)) return -1;
	//    if (!loadRaysFromFile("C:/raydata/barrays128_1.rrfb",rays)) return -1;
	//    if (!loadGeometryFromFile("C:/raydata/frog_neutral.rrf",scene)) return -1;
	//    if (!loadRaysFromFile("C:/raydata/frog_neutralrays128_1.rrfb",rays)) return -1;
	if (!loadGeometryFromFile(sceneFile,scene)) return -1;
	if (!loadRaysFromFile(raysFile,rays)) return -1;

	puts("Building acceleration structure...");
	buildAccelerationStructure(scene);
	Triangle::printSlotStats();

	puts("Shooting rays...");
	float *results = new float[rays.size()];
	shootRays(scene,rays,results);
	shootRays(scene,rays,results);

	const float successRatio = analyzeResults(scene,rays,results,false);
	return (successRatio > 0.995f ? EXIT_SUCCESS : EXIT_FAILURE);
}
